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Abstract 

Bursts of the count rate of extensive air showers (EAS) lead to the appearance 
of clusters in time series that represent EAS arrival times. We apply methods of 
nonlinear time series analysis to 20 EAS cluster events found in the data set obtained 
with the EAS-1000 prototype array. In particular, we use the Grassberger-Procaccia 
algorithm to compute the correlation dimension of the time series in the vicinity of 
the clusters. We find that four cluster events produce signs of chaotic dynamics in 
the corresponding time series. By applying a number of supplementary methods we 
assess that the nature of the observed phenomenon may indeed be chaotic and thus 
deterministic. We suggest a simple qualitative model that might explain an origin 
of EAS clusters in general and "possibly chaotic" clusters in particular. Finally, we 
compare our conclusions with the results of similar investigations performed by the 
E AS-TOP and LA AS groups. 
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1 Introduction 

We have already studied the distribution of arrival times of extensive air show- 
ers (EAS) registered with the EAS-1000 prototype array both by methods of 
classical statistics [1,2] and by methods of cluster analysis [3,4]. In particular, 
we presented 20 EAS cluster events — groups of consecutive showers that were 
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registered in time intervals much shorter than expected ones. These groups 
represent bursts of the EAS count rate and may be considered as cosmic ray 
bursts, see [3,4]. This phenomenon has put forward at least two questions: one 
on the astrophysical nature of the process, another on its statistical proper- 
ties. Namely, while the vast majority of sufficiently long samples satisfy the 
hypothesis for an exponential distribution of time delays between EAS arrival 
times [1,2], the x^-test performed for samples taken in the vicinity of some 
of the EAS clusters made us reject this hypothesis. Thus we have decided to 
apply methods of nonlinear time series analysis to samples that contain EAS 
clusters in order to clarify dynamical reasons of this situation. This approach 
has already proved to be a powerful tool of investigation in different fields of 
science including different fields of astrophysics in general [5] and cosmic ray 
physics in particular, see [6,7,8,9,10,11] and references therein. 

Recall that the data set under consideration represents 203 days of regular 
operation of the array for the period August 30. 1997 to February 1, 1999. 
The total EAS number in the data set equals 1 668 489. The mean number 
of charged particles in a shower is of the order of 1.2 x 10^. This corresponds 
to the energy of a primary particle ~ 1 PeV/nucleon. The mean interval 
between consecutive EAS arrival times is equal to 10.5 s. The discreteness in 
the moments of EAS registration approximately equals 0.055 s (one tic of the 
PC clock). 

In [12], we have already discussed briefly the results of nonhnear analysis of 
one of the EAS clusters. In the present paper, we study this cluster in depth 
and briefly discuss three other clusters that demonstrate nonlinear features. 



2 Preliminary Facts 

In what follows, we study a scalar time series of the form xi,X2, . . . , where 
Xi = ti ~ ti^i, ti is the moment of registration of the ith EAS, i = 1, . . . ,n, 
and to = corresponds to midnight. For the purposes of our investigation, it 
is convenient to use time intervals Xj given in seconds. 

One can point out two main approaches to nonlinear analysis of the time 
series: (1) the scaling exponent and fractal length estimates based on the self- 
similarity properties of experimental signals [13]; (2) an estimate of the corre- 
lation dimension based on the embedding of the original time series in an m- 
dimensional phase space [14,15,16] (see, e.g., [17] for a comprehensive review). 
In the current research, we follow the second approach. Within it, the main 
tool to analyse the dynamics of the time series is the Grassberger-Procaccia 
method [18] with the modification suggested by Theiler [19]. Namely, given 
a sample {xi,X2, ■ ■ ■ ,xn) extracted from the time series we construct m- 
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dimensional delay vectors 

where r is an arbitrary but fixed time increment and m is an embedding 
dimension. Then we compute the number K{p) of vectors with mutual distance 
less or equal than p and such that delay vectors Xj are shifted by at least W 
indices: 

M-W M 

K{p)= E E e(p-||x,-x,||), (1) 

1=1 j=i+W 



where © is the Heaviside step function, M — N — {m— 1)t, and W >lis the 
cut-off parameter (the Theiler window). Finally, we plot log K{p) vs. log p. For 
small p, the slope of this plot is an estimate of the correlation dimension D2: 

dlogp 



where C2 is the correlation sum: 

2K{p) 



C2{p) = 



{M -W){M -W + 1) 



(Obviously, one may use K instead of C2 in the expression for the correlation 
dimension.) A plateau observed in the £)2(p)-plot for small p or, equivalently, 
a so called scaling region in the log K vs. log p plot are regarded as signs of 
chaotic dynamics in the corresponding time series. A value of D2 at the plateau 
is taken as an estimate of the correlation dimension of the attractor underlying 
the data. This quantity also gives a (lower) estimate for the number of degrees 
of freedom in the process under consideration. 

To compute we normally divide each unit interval in Ig p into 50 subinter- 
vals of equal length. We have found that though this number of subintervals 
seems to be small, it adequately reflects the qualitative structure of D2 while 
a much bigger number of subintervals (> 200) leads to considerable fluctua- 
tions in these fluctuations hide the structure of especially for m large 
enough. A smaller number of intervals gives too rough structure of D2- The 
derivative is calculated via a standard three-points algorithm. No smoothing 
or fitting procedures are used. 

Notice that Eq. (1) contains four free parameters: iV, r, W , and m. At the 
preliminary stage of the investigation, we "scanned" the experimental data 
set having split it into adjacent samples with N — 128, 256, 512, and 1024. 
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At this step, we used r = 1, = 1, and odd values of m in the range from 5 
to 13. The value of was chosen to be a power of 2 because this allows one 
to use fast algorithms of calculating the Fourier transform that is the main 
part of traditional Fourier analysis. After we had found a number of samples 
that demonstrated some kind of plateau in the plot of D2{p), we studied the 
corresponding data with different N in the above range and different values 
of r in the range 1-10. (Evidently, in our case r may take only integer values.) 
For A^ large enough, we have performed calculations for m up to 25. In any 
case, m and r were chosen such that the number of delay vectors was greater 
than 100. Finally, to avoid autocorrelation in the time series, we employed 
W^1...20. 

To compute mutual distances between delay vectors (see Eq. (1)), we tried 
several norms: the maximum norm L^, the taxicab norm Li, and the Euclidian 
norm L2. Recall that for a given m-dimensional vector x these norms are 
defined in the following way: 



I'^^lloo l"^*! ' ll-'^lll — l-^il ) I|-''ll2 

l<i<m ,_i 



1=1 



Besides these, we have also employed the "dimension scaled" norms Lie and 
L2C that are expected to provide more reliable results than their classical 
counterparts [20,21]: 




Since the process of calculating D2 is very time consuming, we have used only 
one of these norms, namely the maximum norm, to analyse the complete data 
set. Four other norms were tested for a number of samples both with and 
without a plateau in the D2(p)-plot obtained with Loo. 

There are a number of tools that can help one to verify the results of cal- 
culating the correlation dimension D2, see, e.g., [17]. Among them, we have 
chosen the Theiler-Takens "maximum likelihood" estimator of the correlation 
dimension [22,23]: 

-1 

(2) 




This quantity not only provides an efficient means to find out an optimal 
estimate of the correlation dimension out of the correlation integral but also 
provides an estimate of the statistical error. 
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Another problem in the case when a plateau in a plot of I?2(p) is observed is to 
make an assessment about the nature of the dynamics. The main difficulty is to 
figure out whether one witnesses chaotic dynamics in a deterministic process or 
a special class of stochastic processes [24,25]. One of the main tools to solve this 
problem is the method of surrogate data [26,27], see also [28] and references 
therein. The main idea of this approach is to generate a sufficient number 
of artificial ("surrogate") samples that have the same statistical distribution 
(and possibly some other features) as the experimental data and to study 
the behaviour of the correlation dimension. If surrogate data demonstrate a 
plateau in plots of D2{p) then one concludes that the original data is stochastic; 
otherwise they are likely to be deterministic and chaotic. 

Among other tools one can find the quantity 

C2(m,p) 



suggested in [29] and in the above form — in [30]. As it was shown in [29], 
for a sequence of independent random numbers, C2{m, p) = C2(l,p)™ holds, 
where m is the embedding dimension, and thus iBDs(?7i, p) = 1. Therefore, if 
one finds tsBsinT-, p) 7^ 1 for p where a plateau in the D2(p)-plots is observed 
then one may conclude that the nature of the dynamics is chaotic. 

Besides this, we used a function suggested in [31]: 

= D2(m2,p)-D.(m,,p)^ 

1712 — TTli 



where mi and 1712 are different embedding dimensions. For independent ran- 
dom data, 00 (p) 7^ and 0o — > 1 as p — > 0. We also employed the "normalized 
slope" introduced in [21]: 

(f){m,p) ^ —D2{m,p). (5) 
m 



If 0(m, p) does not converge to in a wide range of p as m grows but to some 
value > 0.1, then most probably data do not represent a chaotic process but 
should be treated by statistical techniques [21]. 

Recall that before applying nonlinear techniques to the analysis of a time se- 
ries, it is strongly suggested to check whether the time series is really nonlinear. 
In particular, one can try a measure for time-reversibility, which is considered 
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to be a good indicator for nonlinearity [32]. For the data sorted in time order, 



is calculated, which is just the mean of the slopes, taken to the third power; 
here Sj are the moments of time and is the variance of the sample. (In our 
case. Si — Si-i — 1 since Xi represent just numerated time delays.) For a time 
series generated by a linear process and for the surrogate data, one expects 
7 fti 0. In contrast, time series with nonlinearities can be asymmetrical in time 
and may yield values of 7 7^ 0. To check, whether 7 significantly deviates 
from zero for the studied sample, one should generate a sufficient number of 
surrogate data. To pay regard to deviations in both directions (7 > and 
7 < 0), a two sided test has to be performed [26]. 

Finally, a few words are in order about stationarity of the time series under 
consideration. As is well known, a fundamental assumption underlying almost 
all existing linear and nonlinear techniques of time series analysis is that the 
time series is stationary, see, e.g., [33]. As we have already mentioned ear- 
lier [2,4,34], the count rate of EAS depends on the value of the atmospheric 
pressure. For the data obtained with the EAS-1000 prototype array, this de- 
pendence can be approximately expressed by a simple formula 

In A^EAS = —PP + const, 

where A^eas is the number of EAS registered in a time unit, P is the atmo- 
spheric pressure, mm Hg, and f3 ~ 10~^ is the barometric coefficient. This 
effect makes the time series non-stationary. To provide stationarity (at large 
time scales) we adjusted time delays between consecutive showers to the at- 
mospheric pressure P* — 742 mm Hg which is close to the average pressure 
for the whole analyzed data set. This adjustment was made by the following 
formula: 



where is the experimental time delay, and Pi is the atmospheric pressure 
at the ith delay. The barometric coefficient was chosen to be /? = 1.08 • 10~^, 
exactly as in our previous paper [4]. Due to the adjustment, the mean and the 
variance of the time series at large time scales were approximately constant. 

Remark 1 In fact, since we shall discuss only comparatively short samples, 
the above adjustment to P* is not important for the present results because at- 
mospheric pressure does not change significantly during periods of time covered 




(6) 



a;, = a;0exp[/5(P*-P,)], 
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Fig. 1. The count rate at a time interval that contains the cluster registered on 
November 11, 1998 (left); the cluster is presented by the central bin. Time de- 
lays (s) between EAS that constitute the cluster (right). The delays are adjusted to 
P* = 742 mm Hg. 

by the samples. We perform the adjustment in order to guarantee stationarity 
of the whole time series. 

To compute the correlation sum C2 effectively, we have worked out an algo- 
rithm based on preliminary sorting of mutual distances between delay vectors. 
Though quite simple, the algorithm occurred to be up to 10 times faster than 
a straightforward computation of C2. To perform calculations, we employed 
GNU Octave [35] running in Mandrake Linux. 



3 The Main Results 



As is well known, one needs sufficiently long samples to perform a successful 
time series analysis. On the other hand, the bursts of the EAS count rate have 
comparatively short time range. Thus we did not in fact expect to find signs of 
chaotic dynamics in the vicinity of EAS clusters. Surprisingly, we have found 
some. 

Having "scanned" the available data set, we found that for the overwhelming 
majority of data samples no scaling region is observed in log C2 vs. log p plots. 
This is quite natural since arrival times of extensive air showers represent a 
simple stochastic process such that for sufficiently long samples, the number 
of EAS registered in a time unit obeys the Poisson distribution in a wide range 
of time delays, see, e.g., [2]. Still, we have come across a number of samples 
with a plateau in the D2(p)-plot in the vicinity of four EAS clusters, namely 
those registered on May 14, November 11, and December 28, 1998 and on 
January 8, 1999. Let us begin our discussion with the second of these events, 
which is the longest one. 

The cluster event observed on November 11, 1998 consists of 136 EAS regis- 
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tered within the period from 01:21:17.47 to 01:38:02.27 (Moscow local time, 
MSK^ ) with the atmospheric pressure P = 752.3 mm Hg, see Fig. 1. Within 
the whole data set obtained on November 11, 1998 the showers that form the 
cluster event have numbers 435-570. The event is made up of three clusters, 
which begin at consecutive arrival times (i.e., the first cluster begins at the 
shower #435, the second one begins at the shower #436, etc.) and end up 
simultaneously at the shower #570. In our opinion, the appearance of three 
clusters does not reflect any process of astrophysical nature but is caused by 
the technique of their selection (see [4] for the details) . Thus we treat the event 
as a single (outer) cluster. The real duration of this cluster equals 1004.8 s 
while the adjusted duration equals 898.8 s; the probability of the appearance 
of such a cluster is of the order of 2 x 10~^. 

Let us take a look at a number of samples in the vicinity of the cluster in order 
to see how the correlation dimension D2 changes when the cluster appears. To 
make this influence more clear we shall not only present samples with a plateau 
in the plot of D2{p) but also a typical sample outside the cluster. Besides 
nonlinear tools, we shall employ the classical Fourier analysis. Recall that it is 
suggested for the Fourier analysis to have x{l) ^ x{N). All samples discussed 
below are chosen to satisfy this demand. In particular, we omit the last shower 
of the cluster and consider a sample that consists of 134 instead of 135 delays. 
Obviously, this does not influence D2 and other quantities signiflcantly. 

Figure 2 shows the correlation dimension (the left column) and the power 
spectrum density (PSD) (the right one) calculated for three samples in the 
vicinity of the cluster. The top row represents a sample that consists of EAS 
#435-569 and thus contains the cluster without the last delay. The length of 
the sample N — 134. The middle row represent a sample that contains the 
cluster but consists of = 256 delays (EAS #425-681). The bottom row 
represents a sample that immediately follows the second one. It also consists 
of 256 delays (EAS #682-938). The correlation dimension was computed for 
r = 1, = 1, and m = 1, . . . , 12. 

What can be seen from this figure? In our opinion, the most important thing 

is that one can see plateaus in both plots of D2{p) that represent samples with 
the cluster. To the contrary, no plateau is observed for the sample outside the 
cluster. We stress that in this sense, the sample is typical for the whole data 
set. 

Next, it is remarkable that the influence of the cluster on the behaviour of 
the correlation is not restricted to the cluster itself, e.g., the sample made of 
EAS #425-681 begins approximately 2.5 min before the cluster and ends up 
in 23.5 min after it. Our analysis has revealed that a plateau can be observed 

^ Recall that MSK is connected to UT as follows: Tmsk = ^UT + 3 + A (hours), 
where A = 1 during Daylight Saving Time, otherwise. 
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Fig. 2. The correlation dimension £>2(p) (the left column) and the power spectrum 
density (the right column) for a number of samples in the vicinity of the cluster 
registered on November 11, 1998. From top to bottom: a sample with the cluster 
(without the last delay), EAS #435-569, N = 134; a sample that contains the 
cluster as a subset, EAS #425-681, N = 256; a sample that follows the second one, 
EAS #682-938, N = 256. The distance p is given in seconds. 

even for samples with ~ 500 lying the in the range of EAS numbers 340- 
1000. It is also interesting to mention that the plateau is more pronounced if 
a sample contains more of the after-cluster showers than those arrived before. 



The level of the plateau changes for different N. The correlation dimension 
fluctuates around D2 ~ 1.77 for the sample with = 134 while D2 ~ 2.45 
for the sample with N — 256. This may be due to low numerousness of the 
samples and/or white noise effects.^ Notice that for both samples with a 
plateau the demand m > 2D2 + 1 is satisfied if we assume that D2 is given 

^ We thank an anonymous referee for pointing this out to us. 
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Fig. 3. The left panel: dependence of the correlation dimension D2{p) on the position 
of the sample {N = 134). The right panel: dependence of the correlation dimension 
at the plateau on the value of the embedding dimension (EAS #425-681, N = 256). 

by the level of the corresponding plateau (see, e.g., [36]). We also remark here 
that the level of a plateau varies for samples with the same length but with 
different positions with respect to the cluster. The left panel in Fig. 3 shows 
the dependence of the mean value of at a plateau on the position of a 
sample. All samples have the same length = 134 but begin at different EAS 
from #425 up to #505. EAS #425 arrived at 01:18:42.04 MSK, EAS #505 
arrived at 01:30:42.15, thus the left ends of these samples cover 12-minute time 
interval. It is evident from the plot that the value of the correlation dimension 
takes different values and monotonically decreases from 1.80 to 1.08. 
For samples that begin at EAS > 510 a plateau becomes shorter and soon 
disappears. 

The right panel in Fig. 3 demonstrates how the value of D2 saturates as 
the embedding dimension m grows. The plot is made for the sample that 
consists of EAS #425-681. It is clearly seen that the saturation takes place 
at m = 9 but not at m = 3 as one may expect for this value of D2 and a 
purely chaotic time scries. As it was shown in [37], this may be due to the 
presence of observational noise and small length of the sample. The demand 
N > 10^2/2 [38] jg aigo obviously hold. 

Finally, one can notice that the plateaus are not as "smooth" and long as 
one may like to see. In this connection, we would like to mention two things: 
(1) Experimental data are affected by noise that can spoil the picture; we plan 
to employ special noise filtering techniques in future. (2) We do not use any 
smoothing techniques while calculating D2. They can improve a plot of D2 
but the situation remains qualitatively the same. 

Now let us turn to the right column of Fig. 2. This column presents the results 
of the Fourier analysis of the three samples. The PSD of the sample without 
a plateau is a broadband one, which is typical for a random noise signal. The 
PSDs of two other samples do not also differ significantly from a broadband 
one. Still it seems worth mentioning that for both samples with a plateau one 
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Fig. 4. The left column: the Theiler-Takens estimator, see Eq. (2). The right column: 
isDS, Eq. (3). The functions are computed for the same samples and values of r, 
W, and m as in Fig. 2. 



can see the highest peaks located at the left end of the spectrum. This is true 
for almost all samples with a plateau in D2-plots we have found. 

Next, let us employ the tools discussed in Section 2 in order to obtain a deeper 
insight in the observed phenomenon. For brevity, let us call the samples with 
a plateau in £)2-plots 'possibly chaotic' (PC). Our aim is to figure out whether 
they are indeed chaotic or represent a special stochastic process. 

As one can see from Fig. 4, the Theiler-Takens estimator (2) does not have a 
clear plateau for either samples but the behaviour of txT for the PC samples 
noticeably differs from that for the third one. Namely, for the PC samples the 
^TT-curves saturate for the corresponding intervals of Ig p for m > 8 while for 



11 



981 1 1 1 : EAS 435...569, AT: 01 :21 :17.47... 01 :37:55.1 8 (N - 134). t - 1 . m - 1 :12, W - 1 . 




Ill: EAS 435...569, AT: 01 :21 :17.47... 01 :37:55.1 8 (N - 134). T - 1 , m - 1 :12. W . 1 . 



981 1 1 1 : EAS 425.. .681 .AT: 01:1 8:42.04, 



1g p 

. 02:01:20.98 (N = 256). T = 1 





0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.5 1.6 

Igp 

981 1 1 1 : EAS 682...938, AT: 02:01 :25.26... 02:46:35.46 (N = 256). T=1,m-1:12.W-1. 




0.8 0.9 1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 



0.8 0.9 1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 



Fig. 5. The left column: the normalized slope (f), see Eq. (5). The right column: the 
function given by Eq. (4), m2 — m\ = 1. The functions are computed for the 
same samples and values of r, W, and m as in Fig. 2. 



the third sample the curves stay separated. 



Next, the right column of Fig. 4 depicts the behaviour of teDSi sec Eq. (3). 
One can see that for both PC samples and m big enough, ^bds noticeably 
differs from 1 at the intervals where plateaus of the correlation dimension are 
observed [Igp G (0.94, 1.12) for the first sample and Igp G (0.94, 1.08) for the 
second one]. Conversely, teos ~ 1 for Igp > 0.94 for the third sample. As 
it follows from the properties of teos discussed above, this kind of behaviour 
witnesses in favour of a chaotic nature of the PC samples and a stochastic 
nature of the sample without a plateau in the D2-plot. 

^ We would like to mention that one can find samples such that sufficiently clear 
plateaus can be observed for both D2- and txT-plots, see the first version of this 
paper for the details [39]. 
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Fig. 5 depicts the normalized slope (see Eq. (5)) and the function (Eq. (4)) 
for the same three samples. One can see that the behaviour of these functions 
for the PC samples noticeably differs from that for the third one. Still, the sit- 
uation is not unequivocal. Namely, for both PC samples, 0o fluctuates around 
zero at the corresponding intervals of Ig p thus giving an argument in favour 
of a chaotic nature of these samples. On the other hand, the limiting values of 
the normalized slope (the bottom curves) lie above 0.1 at the same intervals 
of Ig p thus suggesting a stochastic nature of the samples. This controversy 
can probably be resolved if we recall that converges slowly, see [21]. Really, 
the level of the "hollows" in the plots of decreases if we increase m reaching 
the values ~ 0.1 and ~ 0.15 for the first and the second PC samples respec- 
tively and m = 15. For m > 15 the number K of close delay vectors at the 
corresponding intervals becomes < 100 thus not providing a value statistically 
necessary for accurate computation of the correlation dimension. 

Next, let us say a few words about the results of surrogate data tests. For the 
first PC sample, we generated 99 "shuffled" surrogates obtained by a random 
permutation of time intervals Xi in a given sample. This number of surrogates 
corresponds to a 98% level of significance (L.S.) of the statistical test [32]. For 
the second sample, we made 39 "shuffled" surrogates (95% L.S.) and 39 sur- 
rogates based on the amplitude adjusted Fourier transform method proposed 
in [27]. To obtain these surrogates, we employed the TISEAN package [40]. 
None of these "Fourier-based" surrogates demonstrated a plateau in D2-plots 
thus revealing that the consequence of time delays that constitute the original 
data is crucial for the appearance of plateaus. Therefore, the results of this 
test suggest that the PC samples represent a chaotic (and thus nonlinear) 
process. 

On the other hand, the test for time-reversibility based on 7-statistic (6) did 
not detect nonhnearity in the PC samples. Namely, for the first PC sample 7 = 
0.74 while 7 e (—1.38, 1.28) for the surrogates. For the second PC sample, 7 = 
—0.22 while surrogates gave 7 G (—0.44, 0.45). (To compare, 7 = —0.38 for the 
third sample.) These results imply that the hypothesis for time-reversibility of 
the PC samples cannot be rejected. Therefore, these samples may be linear (or 
obtained by a transformation of a linear process) thus excluding the possibility 
of chaotic dynamics. 

Thus, it is interesting to figure out what kind of stochastic processes can 
describe the distribution of time delays in the PC samples. To do this, we 
performed the x^-test to verify the hypothesis that time delays between con- 
secutive arrival times have an exponential distribution. We used different time 
bins in the range from 1 to 10 s with a step equal to 0.5 s providing that each 
bin taken into account contains more than 10 events. We have found that for 
both PC samples, this hypothesis should be rejected with at least 90% L.S. 
while for the third sample the hypothesis may be accepted with 85% L.S. 
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In addition, we remark that for the above PC samples a plateau was only 
found for r = 1 and the maximum norm. Still there are PC samples in the 
vicinity of the cluster that demonstrate a plateau in the D2-plots for r = 2 and 
both "dimension scaled" norms, see [39] for the details. For all PC samples, a 
plateau is observed in a wide range of the Theiler window W (up to 20). 

Now let us briefly discuss three other EAS clusters that produce signs of 
chaotic behaviour in the corresponding time series. Table 1 contains certain 
information about these clusters as well as the cluster discussed above. Notice 
that while the cluster registered on January 8, 1999 is sufficiently long (A^eas = 
134) two other clusters are short and this makes an appearance of signs of 
chaos in much longer samples very surprising. It is also worth mentioning that 
the fourth cluster (December 28, 1998) was registered within the period of 
observation of the GRB No. 7285 in the BATSE Catalogs [41]. 

Table 1 

EAS clusters that produce signs of chaotic dynamics 



Date 
dd.mm.yy 


Beginning 
hh:mm:ss 


End 
hh:mm:ss 


EAS 
Numbers 


Aeas 




11.11.98 


01:21:17.47 


01:38:02.27 


435-570 


136 


2 X IQ-'^ 


08.01.99 


00:19:47.03 


00:32:26.31 


145-278 


134 


4 X 10"'^ 


14.05.98 


22:24:50.09 


22:24:54.70 


7567-7574 


8 


2 X 10"^ 


28.12.98 


15:22:33.18 


15:23:17.34 


5695-5712 


18 


5 X IQ-'^ 



Notation: Aeas — the number of EAS in a cluster; "Padj — the probability of an ap- 
pearance of the cluster assuming that time delays between consecutive EAS are 
adjusted to the pressure P* = 742 mm Hg and the distribution of the number of 
EAS registered in a time unit (for the whole data set) obeys the Poisson distribution 
with the intensity A ^ 5.77 min"-*^. Moscow local time is given. 

Table 2 presents the results of calculation of the correlation dimension for a 
number of samples in the vicinity of the clusters as well as the results of a 
number of tests discussed above. For the last three samples, the given values 
of D2 were computed for m — 15, 17, and 14 respectively. Notice that the 
PC sample for January 8, 1999 is shifted with respect to the cluster. The 
cluster by itself docs not produce signs of chaotic dynamics. For this event, 
a plateau can be observed not only for r = 1 but also for r = 2 and 3. 
For the cluster registered on May 14, 1998, a plateau can also be observed 
for T = 2. For the last cluster, a plateau was observed only for r = 1. To the 
contrary to most other PC samples, a plateau exists not only for the maximum 
norm but also for the norms L2 and L2c- For all PC samples, a plateau 
was found in a wide range of values of the Theiler window W. This means 
that autocorrelation, which can lead to spurious results, is avoided. 
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Table 2 

Results of different tests applied to PC samples 



Date 


PC EAS 


D2 


max 


^BDS 





00 


SDT 


7-test 


X^-test 


dd.mm.yy 


Numbers 




AD2 














11.11.98 


435-569 


1.77 


0.11 


± 




+ 


+ 








425-681 


2.45 


0.12 






+ 


+ 






08.01.99 


105-232 


0.74 


0.03 


+ 


+ 


+ 


+ 




+ 


14.05.98 


7499-7627 


1.40 


0.04 




+ 


+ 


+ 




+ 


28.12.98 


5632-5760 


1.28 


0.07 


+ 


+ 


+ 


+ 


± 


+ 



Notation: PC EAS Numbers — the range of EAS numbers for a PC sample (to be 
compared with the range of EAS that form the corresponding cluster); D2 — the 
mean value of the correlation dimension at a plateau for a maximum value of m, see 
the text; maxAD2 — the maximum value of the statistical error of D2 {AD2{p) = 
D2{p) / \/ K {p)) . Columns tBDS; 4>i 4>0i SDT (surrogate data test), and 7-test present 
results of the corresponding tests: "+" if a test witnesses in favour of a chaotic 
nature of the sample under consideration, "— " otherwise; "±" means that no definite 
conclusion can be made. The last column presents the results of the x^-test of the 
hypothesis for an exponential distribution of time delays within a sample: "+" means 
that the hypothesis may be accepted for some choice of time bins, "— " otherwise. 
For all data, r = 1, W = 1, and the maximum norm is used. 

As one can see from Table 2, the results of the tests applied to these three 
events are mostly similar to those for the first cluster. Namely, none of the sur- 
rogate samples (39 "shuffled" and 39 Fourier-based surrogates) demonstrates 
a plateau in the D2(p)-plot thus implying a chaotic nature of the original 
samples. On the other hand, the test for time-reversibility usually does not 
detect nonlinearity in the PC samples. The only exception is the last sample 
in Table 2. For this sample, the 7-test applied to shuffled surrogates gives 
the usual result but being applied to Fourier-based surrogates it rejects the 
hypothesis for timc-rcversibility with 95% L.S. It is worth mentioning that 
this result has been confirmed when we employed another technique of mak- 
ing Fourier-based surrogates also implemented in TISEAN. Next, for all PC 
samples the Theiler-Takens estimator demonstrates either a plateau or 
an interval where the curves saturate as the embedding dimension m grows. 
Finally, it is interesting to note that for three PC samples we have found a 
way of grouping time intervals that allows one to accept a hypothesis for an 
exponential distribution of time delays within the corresponding PC sample. 
Still we must remark that for the majority of groupings this hypothesis should 
be rejected with at least 90% L.S. 

Thus we come to a conclusion that the observed dynamics of time intervals 
between consecutive EAS in the vicinity of these four clusters may indeed be 
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Fig. 6. Arrival directions (AD) of EAS that form the cluster registered on Novem- 
ber 11, 1998 and have zenith angle 6" < 40°: •—AD of EAS with AT^ > 9.7 x 10^, 
X — AD of all other EAS. The curve shows the Galactic plane. Equatorial coordi- 
nates are used. 

chaotic though this cannot be assessed unequivocally at the moment. Further 
investigation is necessary. 



4 Discussion 



To get a deeper insight into physics of EAS clusters, let us consider the cluster 
registered on Nov. 11, 1998 from another point of view. A natural question 
that arises when one finds a burst of EAS count rate is whether primary 
particles that led to the appearance of the burst had had a common compact 
source. Fig. 6 shows arrival directions of EAS that belong to the cluster and 
have zenith angle 6 < 40°. Obviously, one cannot give a positive answer to 
the above question. Still, this seems quite natural since the energy of primary 
particles in our case is too low to "remember" the source. Besides this, the 
majority of EAS in the cluster represent a normal count rate. Namely, since 
we register 5.7 EAS/min in average and the duration of the cluster adjusted 
to P* = 742 mm Hg equals approximately 15 min then we may expect to 
register 85-86 EAS within this period. Naturally, arrival directions of these 
EAS should be distributed more or less isotropically. Thus, one may expect 
to observe a kind of a cluster in arrival directions for only remaining 50 EAS. 
Still, this does not happen. One can see doublets and triplets of showers with 
almost coincident arrival directions and a kind of grouping along the Galactic 
plane but no big groups of EAS are observed. 

Another question is whether clusters consist of EAS that are more "powerful" 
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Fig. 7. Model of EAS clusters appearance. Notation: LAP — a layer of accelerated 
particles; ML — magnetic lens. 

than the others. Surprisingly enough, but the answer is negative. As we have 
already mentioned in Section 1, the mean value of electrons iV*"* in EAS in our 
data set is of the order of 1.2 x 10^ particles. To compare, for the above cluster 
TVe ~ 9.7 X 10^ and only 14 EAS have A^e > ^e, see Fig. 6. Two of the_se EAS 
have A^e - 10^ On the other hand, at least 41 EAS have < 1/2N^. (We 
remark that the geometry of the EAS-1000 prototype array does not allow 
one to obtain parameters of all registered EAS.) 

A similar situation is observed for other EAS clusters, both PC and "ordinary" 
ones. Thus, we come to the following observations: (i) EAS in a cluster do not 
have a joint source (arrival direction); (ii) the majority of EAS in a cluster 
have A^'e < iV*"*. This allows us to suggest a simple conceptual model that is 
intended to explain qualitatively how EAS clusters might be produced. To do 
this, we need only to assume that a considerable part of EAS that constitute 
clusters are generated by charged particles. 

As is well known, the interstellar space is filled with strong and highly inhomo- 
geneous magnetic fields. Still, average magnetic fields in the heliosphere are not 
strong enough to influence the dynamics of particles with ii^ > 10 TeV/nucleon 
(see, e.g., [42]). On the other hand, a series of long-term investigations based 
on data obtained with Voyager 1,2, and other space crafts have demonstrated 
that the large-scale magnetic field strength fiuctuations frequently have large 
amplitudes and are intermittent, and that regions of relatively intense mag- 
netic fields can have a radial extent of more than 10 AU, see [43] and references 
therein. 
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Basing on these results, let us consider the following situation. Suppose there 
is an extended and sufficiently "thick" layer of accelerated particles (LAP) 
(possibly called a "wave" ) that moves through space towards Earth, see Fig. 7. 
Normally, if it just passes through, an EAS array registers only a few particles 
originated from the LAP. Now suppose that the LAP meets a region of an 
extended strong and inhomogeneous magnetic field that works as a lens. If 
it happens that this "magnetic lens" declines particles in way that they get 
focused in a "proper" direction then an array may register an excess of EAS 
over the normal count rate, i.e., an EAS cluster. Notice that from an observer's 
point of view showers in the cluster may have very different arrival directions. 
In our opinion, the fact that the majority of EAS in the observed clusters are 
less "powerful" than an average shower witnesses in favour of this model since 
we do not need too strong magnetic fields. The duration of a cluster depends 
on the thickness of the LAP and the time during which the magnetic lens 
"works." 

Now let us ask ourselves how does it happen that the time series that repre- 
sent the majority of EAS clusters (as well as the other data) are stochastic 
but some clusters seem to demonstrate signs of chaotic dynamics. At the mo- 
ment, we cannot give a definite answer but can only conjecture that different 
factors may be involved in this phenomenon. Among them, one can mention 
a possibly fractal nature of the interstellar medium (see, e.g., [44] and ref- 
erences therein), nonlinear mechanisms of particles acceleration, etc. A very 
interesting possibility for EAS time series to become chaotic suggests the fact 
that the large-scale fiuctuations in the interplanetary magnetic field strength 
sometimes have fractal or multifractal structure [45]. This phenomenon has 
been observed at different distances from Earth and for very different time 
scales. It is likely that in case that the "magnetic lens" discussed above has 
certain fractal or multifractal properties then it may lead to the appearance 
of chaotic dynamics in EAS time series. 

Certainly, the presented model is only conceptual and oversimplified. Still we 
hope that it reflects the nature of processes that may lead to the appearance 
of EAS clusters, both "ordinary" and "possibly chaotic" ones. 

Finally, it is interesting to compare our results with the conclusions of similar 
investigations performed by other research groups. In a considerable number 
of articles devoted to the nonlinear time series analysis, one can flnd a com- 
prehensive investigation of EAS arrival times registered with the EAS-TOP 
array [6]. 

Basing on a detailed study of the available experimental data set and the 
results obtained with the underground muon monitor [7] the authors of this 
work made a conclusion that cosmic ray signals are all colour random noise, 
independently of the nature of the secondary particle and of the primary 
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parent particle, but an existence of deterministic chaotic effects in cosmic ray 
time series cannot be completely excluded. It was also demonstrated in one 
of the following articles that an impact of background noise brings additional 
difficulties to the problem of distinguishing between chaotic and stochastic 
dynamics [8]. In our opinion, these conclusions as a whole do not contradict 
our results, especially in view of the fact that signs of chaotic behaviour have 
only been observed for about 0.1% of EAS in the whole data set. 

Besides this, a whole series of investigations devoted to the nonlinear analysis 
of EAS time series are carried out in Japan beginning from early 1990s at 
the experimental arrays that now constitute the LAAS network, see, e.g., 
[10,11] and references therein. The authors of these investigations presented 
several dozens of events that demonstrate chaotic dynamics. More than this, 
it was conjectured that the observed dynamics may be due not only to the 
chaotic structure of the medium through which particles have traversed but 
also to the nature of the primary particles [46]. Later on, there was suggested 
a model according to which chaotic events may be generated by cosmic rays 
that have a structure of a fractal wave arriving from a nonlinear accelerator 
like a supernova remnant [47] . This model needs to be studied in details, but 
seems to be promising. Thus that the results obtained during our analysis 
do not contradict the conclusions of similar investigations performed at other 
EAS arrays. 

The results presented above demonstrate that one can observe an unusual 
dynamics of EAS arrival times in the vicinity of certain clusters of EAS with 
the electron number of the order of 10^. While the overwhelming majority of 
samples in our data set are unambiguously stochastic, a number of samples in 
the vicinity of four EAS clusters demonstrate signs of chaotic dynamics. Still it 
is rather difficult to make a final conclusion on the nature of this phenomenon: 
Does it represent deterministic chaos or a special type of a stochastic process? 
In our opinion, the majority of the tests performed witness in favour of the 
first of these two alternatives. Nevertheless we must mention that our analysis 
may somehow suffer of the fact that the phenomena discussed above are only 
observed at comparatively short time scales with short samples while time 
series analysis usually prefers longer samples. In this connection, we recall 
that our investigation of EAS clusters has revealed an existence of "super 
clusters," i.e., clusters that have duration more than 30 min and consist of 
hundreds of EAS. Our future plans include an analysis of these events. 

It seems to be necessary to continue the work in this area and to involve some 
other methods of nonlinear time series analysis. There are a number of other 
nonlinear tools that may help to make a more definite conclusion about the 
nature of the observed phenomenon. Among them, one can recall space-time- 
separation and recurrence plots and the Lyapunov exponents. There are also 
a number of other measures of nonlinearity besides the one used above [17]. 
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Last but not least, special signal filtering techniques may be used in future to 
reduce effects of background noise on the dynamics of PC samples. 

Finally, as we have already mentioned above, perhaps the biggest puzzle in 
connection with signs of chaos in EAS time series is their astrophysical na- 
ture. It is likely that clusters that produce signs of chaotic dynamics in the 
corresponding time series are similar to the upper part of an iceberg in a sense 
that they do not present the complete process but only the most pronounced 
part of it. We point out that for all PC samples discussed above, the value of 
the correlation dimension is comparatively small {D2 < 3). Since this value 
gives a lower estimate for the number of degrees of freedom in the underlying 
process, it is possible that the structure of this process is not too complicated. 
Still it seems to be a great challenge to work out a good model that could 
explain chaotic dynamics in EAS arrival times. 
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